%% process_boundary_NS_t
function [A, b] = process_boundary_NS_t(A, b, mesh, fespace1, fespace2, pde, t)
Nb1 = fespace1.Nb;
Nlb1 = fespace1.Nlb;
% Neumann
for k = 1:size(pde.bdary,2)
    if 1 == pde.bdary(1,k)
        bde = bdary_edge(mesh, k);
        rows = zeros(size(bde,1)*Nlb1, 1);
        val1 = rows;
        val2 = rows;
        flag = 1;
        for n = 1:size(bde,1)
            points = mesh.node(bde(n,2:3),:);
            AB = [points(2,1)-points(1,1); points(2,2)-points(1,2)];
            vector_n = [0, 1; -1, 0] * AB;
            vector_n = vector_n./norm(vector_n);
            fb1 = @(x,y,t) pde.nu*(pde.u1dx(x,y,t)-pde.p(x,y,t))*vector_n(1) + pde.nu*pde.u1dy(x,y,t)*vector_n(2);
            fb2 = @(x,y,t) pde.nu*pde.u2dx(x,y,t)*vector_n(1) + pde.nu*(pde.u2dy(x,y,t)-pde.p(x,y,t))*vector_n(2);
            e = bde(n,1);
            V = mesh.node(mesh.elem(e,:),:);
            [int_points_loc, int_weights_loc] = integral_loc(points, "1d");
            phi = basis_loc(V, int_points_loc(:,1), int_points_loc(:,2), fespace1.basis_type, "x");
            ffb1 = fb1(int_points_loc(:,1), int_points_loc(:,2), t);
            ffb2 = fb2(int_points_loc(:,1), int_points_loc(:,2), t);
            for i = 1:Nlb1
                f1 = ffb1 .* phi(:,i);
                f2 = ffb2 .* phi(:,i);
                r1 = int_weights_loc * f1;
                r2 = int_weights_loc * f2;
                rows(flag) = fespace1.elem(e,i);
                val1(flag) = r1;
                val2(flag) = r2;
                flag = flag + 1;
            end
        end
        b1 = sparse(rows, 1, val1, Nb1, 1);
        b2 = sparse(rows, 1, val2, Nb1, 1);
        vector_b = [b1; b2];
        b(1:2*Nb1) = b(1:2*Nb1) + vector_b;
    end
end
% Dirichlet
for k = 1:size(pde.bdary,2)
    if 0 == pde.bdary(1,k)
        bdn1 = bdary_node(mesh, fespace1, k);
        A(bdn1,:) = 0;
        A(sub2ind(size(A), bdn1, bdn1)) = 1;
        b(bdn1) = pde.u1(fespace1.node(bdn1,1),fespace1.node(bdn1,2),t);
        bdn2 = bdn1 + Nb1;
        A(bdn2,:) = 0;
        A(sub2ind(size(A), bdn2, bdn2)) = 1;
        b(bdn2) = pde.u2(fespace1.node(bdn1,1),fespace1.node(bdn1,2),t);
    end
end
% stablization
% M = term_P_Q(mesh, fespace2);
% matrix_A = blkdiag(sparse(2*Nb1,2*Nb1), 1e-6*M);
% A = A + matrix_A;
A(end, :) = 0;
A(end, end) = 1;
b(end) = pde.p(fespace2.node(end,1),fespace2.node(end,2),t);
end
